1987 


N88- 1560 


NASA/ASEE SUMMER FACULTY RESEARCH FELLOWSHIP PROGRAM 


MARSHALL SPACE FLIGHT CENTER 
THE UNIVERSITY OF ALABAMA IN HUNTSVILLE 


INVESTIGATION OF THE FEASIBILITY OF AN ANALYTICAL 
METHOD OF ACCOUNTING FOR THE EFFECTS OF ATMOSPHERIC 
DRAG ON SATELLITE MOTION 


Prepared By: 

Robert E. Bozeman 

Academic Rank: 

Professor 

University and Department 

Morehouse College 
Mathematics Department 

NASA/MS FC 

Systems Analysis and 
Integration 

Laboratory : 

Division : 

Mission Analysis 

Branch : 

Orbital Mechanics 

NASA Colleague: 

Larry D. Mullins 

Date : 

August 7,1987 

Contract No: 

The University Of 
Alabama in Huntsville 
NGT-01-008-021 




Vll-i 



ABSTRACT 


An analytic technique for accounting for the joint 
effects of earth oblateness and atmospheric drag on 
close-earth satellites is investigated. The techni- 
que is analytic in the sense that explicit solutions 
to the Lagrange planetary equations are given; con- 
sequently, no numerical integrations are required in 
the solution process. The atmospheric density in the 
technique described in this report is represented by 
a rotating spherical exponential model with superposed 
effects of the oblate atmosphere and the diurnal vari- 
tions. A computer program implementing the process is 
discussed and sample output is compared with output 
from program NSEP (Numerical Satellite Ephemeris Pro- 
gram). NSEP uses a numerical integration technique to 
account for atmospheric drag effects. 
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INTRODUCTION 


In order to accurately predict orbital lifetimes and obtain 
trajectory data on close-earth satellites, it is generally 
accepted that joint effects due to pertubations such as 
atmospheric drag and the oblateness of the earth must be 
included in the mathematical model. However, inclusion of 
these joint pertubative effects usually leads to coupled 
systems of nonlinear differential equations that are diffi- 
cult to solve in closed form - i.e., obtain an exact solu- 
tion. Consequently, most computational procedures that 
attempt to account for joint effects make use of numerical 
integration techniques to generate solutions to the coupled 
equations. Although these procedures can produce accurate 
results, they are often found to be time-consuming in their 
execution because of the numerical integrations involved. 
Because of this, persons working in real-time satellite 
tracking activities are continually seeking better ways of 
accelerating the solution process. 


This report is an attempt to evaluate and test an analytic 
model that includes the joint effects of oblateness and 
drag and requires no numerical integration in its execu- 
tion. The theory for this model was developed by Sean C. 

H. Chen ([4], [5]). Chen's development commences with the 
Lagrange Planetary equations for the classical orbital ele- 
ments. From the outset the development considers jointly 
the coupled effects of oblateness and drag and is not the 
superposition of two separate developments. Consequently, 
it is considered to be a mathematically rigorous theory. 

The solutions to the variational equations of motion are 
derived using a two-variable asymptotic expansion. Per- 
tubations due to the second harmonic are considered as 
first order effects and pertubations due to the third and 
fourth harmonics and the drag force are considered as 
second order effects for the series expansion . The atmos- 
pheric density model used in the development of the solu- 
tions is a rotating spherical exponential model with super- 
posed effects of the oblate atmosphere and the diurnal 
variation . 


A computer program based on the drag theory of Chen is 
described and sample output is discussed. 
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GENERAL FORM OF SOLUTIONS 


The general form of the solutions to the equations of 
motion for the drag problem will be briefly described. 
The solutions themselves are quite lengthy and will not 
be included in this document. The interested reader can 
consult references [4] and [5] for complete details. 


The equations of motion of a particle of mass m, under 
the attraction of a mass m ,but acted on by perturbing 
forces F = (R,S,W), can be derived in terms of the classi- 
cal orbital elements - a , e , i ,Cl , M. Here, a is the semi- 
major axis; e is the eccentricity? i is the inclination of 
the orbit; 1 ~L is the right ascension of the ascending node; 
( 4 ) is the argument of perigee; and M is the mean anomaly. 
The Lagrange planetary equations for the orbital elements 
can be expressed in the general form: 


da/dt 
de/dt 
di/dt 
d -^ 2 /dt 
d^/dt 


2 

= ((2a )4^p) [ e ( sinv ) R + (p/r)S] 

= I (sinv)R + ({( r+p ) cosv + er}/p)S] 

= ( ( rcosu ) /^U P ]W (1) 

= [ ( rsinu sini)]W 

= ( 1/e ) (/p/& )[-( cosv ) R + ( r/p + l)(sinv)S - 

( ( re( sinu) ( coti ) )/p}W] 


dM/dt = n- (1/^Ta )[ (2r-(p/e)cosv)R + { ( ( r+p ) s inv ) /p } S 

In the above expressions, r is the magnitude of the 
position vector between the two bodies, n is the mean 
motion, v is the true anomaly, u is the argument of latitude, 
p = a ( 1 _ e* ), and A-= k(m+m ), (k = gravitational constant). 


In deriving solutions to the planetary equations, Chen 
assumes that the perturbing forces F are of the form 
p = p ( c ) + F ( D ) where F(C) are the conservative forces 
which includes the earth's oblateness and F(D) is the 
drag force. The drag force per unit mass is represented by 

F ( D ) = -0.5CD(A/m)j) VV (2) 

where CD is a nondimensj^onal drag coefficient, m is the 
mass of the satellite, V is the velocity of the satellite 
relative to the atmosphere, A is the effective cross 
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sectional area of the satellite and P represents the local 
atmospheric density. The analytic model for the density is 
taken to be: 




exp( - ( h-h 0 )/H ) [ 1 + F ( Acosv + 
F = ( RHOl— RH02 )/( RH01+RH02 ) 

= RHOl/( 1 + F ) 

H = -2ae/ln( RH03/RH01 ) 


Bsinv ) ] 


(3) 


where RHOl is a value of the density computed at the center 
of the diurnal bulge, RHOl is a density value computed at 
perigee hight at a point opposite the bulge and RH03 is a 
density value computed at apogee height at a point under 
the bulge. A,B,h, and h d are parameters that depend on the 
position of the point under consideration. 


The Lagrange planetary equations (system (1)) with F(C) 
replaced by the earth gravitational potential up to and 
including the fourth harmonic and with F(D) replaced by 
(2) and (3) are solved analytically for the orbital 
elements. A two variable asymptotic expansion technique 
is used to generate the solutions. The solution for each 
orbital element is obtained and expressed in the form: 

X = X( OBLATE) + X ( DRAG ) 

where X represents any of the orbital elements, X ( OBLATE ) 
is the part of the solution generated by the conservative 
forces and X ( DRAG ) is the portion generated by the drag 
effects. The solutions are quite lengthy and will not be 
reproduced here. The X(OBLATE) components are contained 
in reference [5] (pages 6-3 -6-7) and the X ( DRAG ) compo- 
nents are contained in reference [4] (pages 6-1 - 6-7) . 
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PROGRAM DESIGN AND SAMPLE OUTPUT 


The Analytic Satellite Ephemeris Program (ASEP) now in use 
in the Orbital Mechanics branch of the Mission Analysis 
Division at Marshall Space Flight Center is designed to 
accurately simulate the motion of a satellite in earth 
orbit (see [10]). The solution to the satellite motion is 
completely analytic and requires no numerical integrations 
for its execution. The program is based on the theoretical 
work of Sean C.H. Chen as described in [5] and only in- 
cludes pertubations due to the oblate earth. Atmospheric 
drag effects are not included. In designing a test pro- 
gram that included both oblateness and drag pertubations, 
a subroutine (ORBIT) from the ASEP program was selected for 
modification. Subroutine ORBIT accepts initial values for 
the orbital elements as input and generates a detailed 
trajectory. In designing the program, several tasks were 
executed : 

(1) The mathematical development of the analy- 
tic solution derived by Chen and con- 
tained in [4] and [5] was reviewed; 

(2) All subroutines in ORBIT that were impacted 
by the drag terms were identified and 
appropriately modified; 

(3) A driver (main calling program) for ORBIT 
was written. 

As a result of the above tasks, a test program resulted. 

The components of the program and a brief description of 
each are listed below: 

BDRAG - the main driver program. 

BORBIT - BORBIT is the ORBIT subroutine from 
ASEP containing minor changes to re- 
flect the inclusion of the drag terms. 

BRATES - BRATES contains modifications to sub- 
routine RATES of ASEP that are used to 
generate the secular terms in the 
solution of the orbital element. 

BSECTM - BSECTM is used together with BRATES to 
generate the secular terms. 

BLPTRM - BLPTRM contains modifications to sub- 
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routine LPTERM of ASEP that are used 
to generate the long period components 
of the solution. 

FIT - Program FIT evaluates the constants H, 

F and , given in (3). These values 
are needed for the local density. 


The drag terms are primarily accounted for in the sub- 
routines BRATES, BSECTM, and BLPTRM. The modifications 
are quite extensive and are included in the Appendix. The 
interested reader can compare the expressions in the three 
subroutines with the solutions for the orbital elements 
contained in [4] and [5]. 


The test program described above has been coded and im- 
plemented and the results compared with output from NSEP 
(Numerical Satellite Ephemeris Program). NSEP handles 
the drag effects by a numerical integration process. The 
NSEP program has been used extensively and found to be 
accurate. However, the implementation is slow because of 
the numerical integrations involved. 


Results from the analytic model (no numerical integra- 
tions required) described in this document are compared 
with NSEP output. A comparison of agreement between 
the two procedures is made simply by subtracting the 
corresponding values given by the two techniques for 
the orbital elements. A difference of zero would in- 
dicate that both procedures gave the same value for that 
element at the given time. Charts #1 and #2 give the de- 
viations in five orbital elements over a ten day period. 
It can be noted from the charts that there is reasonably 
good agreement for the eccentricity, argument of perigee, 
and right ascension of the ascending node. However, both 
samples indicate less than favorable comparisons in the 
semi-major axis and the mean anomaly. It is believed 
that better comparisons will result once a good fit of 
the parameters H, F, and 0 needed for the atmospheric 
density, is obtained. * c 
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COMPARISON WITH NSEP OUTPUT 
[ X ( NSEP ) -X ( ANALYTI C ) , X= ORBITAL ELEMENT] 


p = 5 . 80595D-11 
'Q 

INITIAL CONDITIONS: 


INPUT DATA 
F = -.01599 


a = 6755.734(km) 
O. = 3 . 77821 ( deg ) 
i = 57 . 0 ( deg ) 


H = 49797. 8(M) 


e = .019193 
CO = 180 . 0 ( deg ) 
M = 0 . 00 ( deg ) 


ORBITAL ELEMENT 


MISSION ELAPSED TIME (DAYS) 


4 

6 

8 

10 


SEMI-MAJOR AXIS(KM) 

-.575 

-.868 

-1 . 168 

-1 . 467 

ECCENTRICITY 

-.00032 

-.00046 

-.00063 

-.00076 

RIGHT ASCENSION OF 
ASCENDING NODE ( DEG ) 

.03 

.04 

. 06 

.07 

ARGUMENT OF PERIGEE 
(DEG) 

.03 

.03 

.17 

.28 

MEAN ANOMALY (DEG) 

4 . 50 

10.13 

17.89 

27.93 


CHART #1 


COMPARISON WITH NSEP OUTPUT 


[ X ( NSEP ) — X ( ANALYTIC ) , X= ORBITAL 

ELEMENT] 


INPUT DATA 


D = 4 . 45D-11 F = 

0.99 H = 1039000 ( M) 

Jo 



INITIAL CONDITIONS: a 

= 6755.734(km) 

e - .019193 

.12 

= 3 . 77821 ( deg ) 

cO = 180 . 0 ( deg 

i 

= 57 . 0 ( deg ) 

M = 0.00(deg) 


ORBITAL ELEMENT 


MISSION 

ELAPSED TIME 

( DAYS ) 


4 

6 

8 

10 


SEMI-MAJOR AXIS (KM) 

1.635 

2.480 

3 . 326 

4.182 

ECCENTRICITY 

-.00024 

-.00034 

-.00046 

-.00053 

RIGHT ASCENSION OF 
ASCENDING NODE ( DEG ) 

.04 

. 06 

. 09 

.12 

ARGUMENT OF PERIGEE 
( DEG) 

-.27 

-.39 

-.41 

-.43 

MEAN ANOMALY (DEG) 

4.75 

10.47 

18 .33 

28.45 


CHART #2 


I 
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CONCLUSIONS 


An analytic model for accounting for joint effects of 
earth oblateness and atmospheric drag on close-earth 
satellites has been described. Preliminary tests in- 
dicate reasonable results when compared with output 
generated by the NSEP program (Numerical integration 
routine). However, more test cases will be required 
before a final determination can be made on the accu- 
racy of the analytic process. Also, it will be useful 
to compare output from the analytic method with actual 
data collected from observations of close-earth satel- 
lite motion. Such data was not available during the 
period of this report. 


A number of authors have pointed out the difficulty 
in describing a realistic model of the atmospheric 
density because it depends on position and time in a 
very complex manner. In the model tested in this 
report, it was necessary to fit three constants F, / 0 , 
and H in order to describe the atmospheric density. 
This fitting was done by first using the Jacchia model 
to calculate values of the density at specific points 
in the diurnal bulge. These values are used to com- 


pute the numerical values of the three parameters. 
More research is required on this process so that a 
determination can be made on when a " good " fit 


has been obtained. Additionally, in the tests conduct 


ed in this study, only one fitting of the constants 
was done for the entire ten-day mission. It is an- 
ticipated that better results will be obtained if the 
constants were updated at least every twenty-four 
hours. This updating process needs to be tested. 
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c BR A T E S 



c 

c This routine computes the following quantities : 

c 

c 1. Mean Motion 

c 2. Perturbed Mean Motion 

c 3. Nodal regression rate 

c 4. Rate of rotation of the Line of Apsides, 

c 5. Lambda dot (one of the regularized elements) 

c 6. Nodal Period 

c 7. Kepplerian Period 

c 

+ 

c 


SUBROUTINE BRATES 

DOUBLE PRECISION PI,TPI,P02,0ME,AMU,RE,P0LR,AMI,P0L,0E,ME, 

*IME,CAR,EMI,APMI,RANMI,MAMI,UMI,TAMI,LAMMI,IMI,PM,AK2,AAK2, 

*AK4, PERK, PERN, MMOT,LAMDOT,PMMOT,APDOT,RANDOT, 

*CD,A,DM,RP,B,H,FE,QO,DRG,RHONUL,G,DB,DBSQR,DBCUBE,ESQR, 

*DD,EL,Q,B0,B1,YA,ZA,YE,ZE,YI,ZI,PA,PE,DPI,P1A,P1E,YISTAR, 

*ZISTAR,PISTAR,Y1ISTR,ZIISTR,Y3I,Z3I,Q1ISTR,Q3I,P1I,TDG, 

*TPIOT, ELM, BB2,DADOT,DEDOT, DIDOT, DRADOT,DAPDOT,DMADOT, 
*MBF0,MBF1 ,F, ALI 
REAL NODPER,KEPPER 

C0MM0N/BL0K1/PI,TPI,AMU,RE,P0LR,AJ2,AJ3,0ME,P02,AJ4 
COMMON/BLOK8/AMI,EMI,IMI,APMI,RANMI,MAMI,UMI,TAMI,LAMMI 
COMMON/BLOKll/MMOT,PMMOT,APDOT,RANDOT,LAMDOT 
COMMON/BLOK 1 2/JCIRC 

COMMON/BLOK 1 7/POL(6) , OE(6) , ME(6) , IME(6) , CAR(6) 
COMMON/BLOK57/NODPER,KEPPER 

COMMON/BLOKllO/DADOT, DEDOT, DIDOT, DRADOT,DAPDOT,DMADOT 

COMMON/BLOK1 1 1/PA,PE,DPI,PISTAR,Q3I 

COMMON/BLOKI13/DRG,F,RHONUL,H,ALI 

DATA DPR/57. 29577951D + 0/RPD/.0174532925D + 0/EL/.0000729D + 0/ 

DATA B/6.356780D + 6/ 

CALL UNSIME 
OME2 = 1 .-EMI**2 
BB = .5*SIN(IMI)**2 
AA= 1./3.-BB 
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PM = AMI*OME2 

MMOT = DSQRT(AMU/AMI* *3) ! mean motion 

AK2 = 1 . 5*AJ2*(RE/PM)**2 
AAK2 = AK2*(PM/AMI)**2 
AK4 = AJ4/(AJ2**2) 

RANDOT = -AK2*MMOT*COS(IMI)*(1.-AK2*(-1.5+10./3.*BB + EMI**2 
**(-l./6.-5./12.*BB)-3.*AA*SQRT(OME2) + 35./18.*AK4*(-6./7. +3.*BB) 

**(1. + 1.5*EMI**2))) ! Nodal Regression rate 

APDOT = AK2*MMOT*((2.-5.*BB) + AK2*(4.-103./6.*BB + 215./12.*BB**2 + 
*EMI**2*(7./12.-.75*BB-15./8.*BB**2) + (AA+ 15.*AA**2)*SQRT(OME2) 

* + 35./18.*AK4*(12./7.-93./7.*BB + 21.*BB**2 + EMI**2*(27./14.-189./ 

*14.*BB + 81./4.*BB**2)))) ! Rot. rate of line of apsides 

PMMOT = MMOT* ( 1. +AK2*SQRT(OME2)*(( 1.-3. *BB) + AK2*(1 1./6. -26. /3.*BB 

* + 125./12.*BB**2 + EMI**2*(2./3.-4./3.*BB-5./3.*BB**2) + 21./2.*AA**2 
**SQRT(OME2) + 35./18.*AK4*EMI**2*(9./14.-45./7.*BB + 45./4.*BB**2))) 

*) ! perturbed mean motion 

LAMDOT = MMOT*(l.+AAK2*((3.-8.*BB)-AAK2*(35./6.-155./6.*BB + 85./3.* 
*BB**2 + AA + 51./2.*AA**2 + 35./18.*AK4*(12./7.-93./7.*BB + 21.*BB**2))) 

*) ! Lamda dot (regularized element) 

PERN = (TPI/(PMMOT + APDOT))/3600. ! Nodal Period 


NODPER = PERN 

PERK = TPI/MMOT/3600. ! Kepplerian (2-body) Period 

KEPPER = PERK 
C 

D TYPE VRHONUL = ’,RHONUL 

D TYPE *,’F = \F 

D TYPE *,’H = \H 

RP = AMI*(1.-EMI) 

FE = (RE -B)/RE 

QO = RE*FE*SIN(IMI)* :,: 2*SIN(APMI)**2/H 
G = DRG*RHONUL*EXP(-(AMI-RP)/H +Q0) 

DB = AMI*EMI/H 
D TYPE VDB = ’,DB 
DBSQR = DB**2 
DBCUBE = DBSQR*DB 
ESQR = EMI* *2 

DD = EL/MMOT*SQRT(OME2)*COS(IMI) 

Q = (RE-B)/H*SIN(IMI)**2 

BO = MBFO(DB) 

B1 = MBFl(DB) 
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DEFINITION OF PK (K = A,E,I) 


YA = (l.-DD)**2 + EMI**2*(1.5 + DD) 

ZA = 2.*EMI*(1 .-DD**2)-ESQR/DB*(1 .5 + DD) 

YE = EMI*(1.-DD)*(1. +3.*DD)-EMI**2/(2.*DB) 

ZE = (l.-DD)**2 + EMI**2*(0.5-DD)-EMI/(2*DB)* 

1 (l.-DD)*(2. + 5.*DD) + EMI**2/DB**2 

YI = (l.-DD) + EMI**2/2.*(l.-13.*DD) 

ZI = -2.*EMI*(1.-2.*DD)-EMI**2/(2.*DB)*(1.-13.*DD) 

PA = YA*B0 + ZA*B1 

PE = YE*B0 + ZE*B1 

DPI = YI*B0 + ZI*B1 

DEFINITION OF P1A,P1E,AND P1I 

P1A = -0.5*(l.-3./8.*Q)*PA 

PIE = -0.5*(l.-3./8.*Q)*PE 

YISTAR = (1. -DD) + ESQR/2. *(1.-13. *DD)-4.*EMI*DD/DB- 
1 3.*ESQR*(1.+3.*DD)/DBSQR 

ZISTAR = -2.*EMI*(1.-2.*DD)-1./DB*(2.*(1.-DD)-2.5*ESQR* 

1 (1. +3.*DD)) + 8.*EMI*DD/DBSQR+6.*ESQR*(1. +3.*DD) 

2 /DBCUBE 

PISTAR = YISTAR*B0 + ZISTAR*B1 

Y1ISTR = l./DB*(2.*EMI*DD-6./DB*((l.-DD)-.75*ESQR*(5.-DD)) + 
1 48.*EMI/DBSQR) 

Z1ISTR = 1 ./DB*(1 .-DD)-ESQR/2.*(1 1 . + DD)-4. *EMI*(3. + DD)/DB + 
1 12./DBSQR*((1.-DD)-.75*ESQR*(5.-DD))-96.*EMI/DBCUBE 

Y3I = l./DB*(2.*EMI !,! DD-3./DB*(l.-DD)^(l. + ESQR) + 

1 24. *EMI/DBSQR-60. *ESQR/DBCUBE*(2. 5-. 5*DD)) 

Z3I = l./DB*((l.-DD)-ESQR/2.*(l.+3.*DD)-2.*EMI/DB* 

1 (3. +2.*DD)+3./DBSQR*(2.*(l.-DD) + ESQR/2.*(29.-9.*DD)) 

2 -48.*EMI/DBCUBE + 60.*ESQR*(5.-DD)/DB**4) 

Q1ISTR = (1.-ESQR)*(Y1ISTR*B0 + Z1ISTR*B1) 

Q3I = (l.-ESQR)*(Y3I*B0 + Z3I*Bl) 

P1I = -0.5*(1.-.5*Q)*(DPI-.5*PISTAR+Q1ISTR + 2.*Q3I) 

COMPUTATION OF ORBITAL ELEMENTS WITH DRAG 


TDG = TPI/MMOT*(l.-AK2*SQRT(OME2)*(l.-3.*BB)) 
TPIOT = TPI/TDG 
ELM = EL/MMOT 
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BB2 = SIN(IMI)**2 

DADOT = -TPIOT*(2.*G*AMI**2*(PA + Q*PlA)) 

DEDOT = -TPIOT*(2.*G*AMI*OME2*(PE + Q*PlE)) 

DIDOT = -TPIOT*(. 5*G*AMI*ELM*SIN(IMI)/SQRT(OME2)*(DPI + Q*P I 1)) 
DRADOT = -TPIOT*(G*AK2*MMOT*AMI*(COS(IMI)*((7./2.*PA-4.*EMI* 

1 PE) + Q*(7./2.*PIA-4.*EMI*P1E)) + .25*ELM*BB2/ 

2 SQRT(OME2)*(DPI + Q*PlI))) 

DAPDOT = TPIOT*(G*AK2*MMOT*AMI*((2.-5./2.*BB2)*((7./2.*PA-4.* 

1 EMI*PE) + Q*(7./2.*PlA-4.*EMI*PlE)) + 5./4.*ELM* 

2 BB2*COS(IMI)*(DPI + Q*P1 1))) 

DMADOT = TPIOT*(G*AK2*MMOT*AMI*(SQRT(OME2)*(l.-3.*BB)* 

1 ((7./2.*PA-3.*EMI*PE) + Q*(7./2.*PlA-3.*EMI*PlE)) 

2 + .75*ELM*BB2*COS(IMI)*(DPI + Q*P1I))) 

RETURN 

END 
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c- — 
c 

r 

BSECTM 

+ 

V/ 

c 

c 

This routine computes the Secular 

terms of the solution 

c 

to the satellite motion. 


c 

c 


b 


SUBROUTINE BSECTM 

DOUBLE PRECISION T,TO,TF,TAPO,TPER,GET,AMI,EMI,IMI,APMI, 
*MAMI,UMI,TAMI,LAMMI,GSTC,JULDAY, APSEC, RANSEC, MASEC, 
*MMOT,PMMOT,APDOT,RANDOT,LAMDOT,PRIN,DADOT, DEDOT, 
*DAPDOT,DMADOT, DRASEC, DAPSEC, DMASEC, DASEC, DESEC, DISEC, 
*RANMI,LAMSEC, DIDOT, DRADO 

COMMON/BLOK2/T,TO,TF,TAPO,TPER,GET,DAYNO,JULDAY,GSTC 

COMMON/BLOK8/AMI,EMI,IMI,APMI,RANMI,MAMI,UMI,TAMI,LAMMI 

COMMON/BLOKll/MMOT,PMMOT,APDOT,RANDOT,LAMDOT 

COMMON/BLOK12/JCIRC 

COMMON/BLOK20/APSEC,RANSEC,MASEC,LAMSEC 
COMMON/BLOK11 0/D ADOT, DEDOT, DIDOT, DRADOT,DAPDOT,DMADOT 
COMMON/BLOK112/DASEC, DESEC, DISEC 
CALL UNSIME 

RANSEC = RANMI + RANDOT* (T-TO) 

DRASEC = DRADOT*(T-TO)**2 
RANSEC = PRIN(RANSEC + DRASEC) 

IF(JCIRC.EQ.2) GO TO 1 

LAMSEC = PRIN(LAMMI + LAMDOT* (T -TO)) 

RETURN 

1 APSEC = APMI + APDOT* (T-TO) 

DAPSEC = DAPDOT*(T-TO)**2 
MASEC = M AMI + PM MOT* (T-TO) 

DMASEC = DMADOT*(T-TO)**2 
DASEC = D ADOT* (T-TO) 

DESEC = DEDOT* (T-TO) 

DISEC = PRIN(DIDOT* (T -TO)) 

APSEC = PRIN(APSEC + DAPSEC) 

MASEC = PRIN(MASEC + DMASEC) 

RETURN 

END 
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c 
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+ 

B L P T R M 

+ 

This routine computes the Long Period Terms that are a 
part of the satellite ephemeris solution equations. 

+ 


SUBROUTINE BLPTRM 

DOUBLE PRECISION 0ME,P02,PI,TPI,AMU,RE,P0LR,T,T0,TF,TAP0, 
*get,pol,oe,me,ime,ami,pm,reoa,reop,emi,imi,apmi,ranmi, 
*UMI,TAMI,LAMMI,BB,ALP,BO,ABET,QOB,EOB,APSTAR,MMSTAR, 
*ELP,RANLP,MALP,LAMLP,ETAB,ZMB,AK42,AK2,AK32,MMOT, 
*APDOT,RANDOT,LAMDOT,ETAI,TPER,MAMI,ILP,APLP,PMMOT, 
*CD,A,DM,RP,B,H,FE,QO,DRG,RHONUL,G,DB,DBSQR,DBCUBE, 
*Q,B0,B1»MBF0,MBF1 ,PA,PE,DPI,PISTAR,Q3I,TDG,TPIOT,ELM,EP, 
*EB,DB1,EBSQR,Y1A,Z1A,Y1E,Z1E,Y1I,Z1I,Y1ISTR,Z1ISTR,Q1A,01E, 
*Q1I,Q1ISTR,P2A,P2E,P2I,P3I,Y3W,Z3W,Y3M,Z3M,Q3W,Q3M,P2W, 
*Y5A,Z5A,Y5E,Z5E,Y5I,Z5I,Y5ISTR,Z5ISTR,P5A,P5E,P5I,P5ISTR,Y6I, 
*V6I,Y6W,U6W,Z6W,V6W,U6M,Z6M,V6M,P6I,P6W,P6M,EPSl,ALDOT, 
*SEPS1 ,COISQ,SISQ,APDT,RADT,ELl,EL2,EL3,EL4,EL5,EL6,EL7,EL8, 
*SLW,SLWO,DLW,DLWO,SLOW,SLOWO,DLOW,DLOWO,SL3W,SL3WO, 
*SLO3W,SLO3W0,DLO3W,DLO3W0,SNA,SNA0,SNB,SNB0,SNC, 
*SNCO,SND,SNDO,SNE,SNEO,SNF,SNFO,SNG,SNGO,SNH,SNHO,CSA, 

*csb,csbo,csc,csco,csd,csdo,awbar,aprim,bprim,acprm, 

*ADLPM,AOPRIM,BOPRIM,ACOPRM,BSOPRM,AODLPM,DAPM,DS2W, 

*DS4WI,DC2W,DC2WI,DSI,DSISQR,GOEP,REOPSQ,ALPD,ELPD,DILP, 
*APLPD,DMALP,U6I,Z6I,F,ALI,ESQR,DD,EL,P2M,DLMDA,DL3W, 
*DL3W0,CSA0,BSPRM,DS2WI,DS4W,RANLPD 
COMMON/BLOK 1 2/JCIRC 

C0MM0N/BL0K1/PI,TPI,AMU,RE,P0LR,AJ2,AJ3,0ME,P02,AJ4 

COMMON/BLOK2/T,TO,TF,TAPO,TPER,GET,DAYNO,JULDAY,GSTC 

COMMON/BLOK17/POL(6),OE(6),ME(6),IME(6) 

COMMON/BLOK14/BB,ALP,BO,ABET,QOB,EOB 

common/blokii/mmot,pmmot,apdot,randot,lamdot 

COMMON/BLOK21/ELP,ILP,APLP,RANLP,MALP,LAMLP,ALPD 
COMMON/BLOK8/AMLEMUMLAPMLRANMLMAMLUMLTAMLLAMMI 
COMMON/BLOK11 1/PA,PE,DPI,PISTAR,Q3I 
COMMON/BLOK113/DRG,F,RHONUL,H,ALI 
DATA PI/3. 141 592654D + O/TPI/6.283 185307D + 0/ 
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DATA AJ2/1. 0827E-3/AJ3/-2. 56E-6/AJ4/- 1 . 58E-6/ 

DATA RE/6378160.D + 0/AMU/3.986012D + 14/EL/. 0000729D + 0/ 

DATA DPR/57. 29577951D + 0/B/6.356780D + 6/RPD/.0174532925D + 0/ 

CALL UNSIME 

AA= l./3.-.5*SIN(IMI)**2 

BB= 1./3.-AA 

OME2 = 1 .-EMI*EMI 

PM = AMI*OME2 

C TYPE 2002,APMI,AJ2,RE,EMI,BB,MMOT,T,TO 

C2002 FORMAT(D20. 12,E20.7,/,3D20. 12,3D20. 12) 

APSTAR = APMI + 1 . 5*AJ2*(RE/PM)**2*(2.-5. * BB)*MMOT* (T-TO) 
MMSTAR = MAMI + MMOT*(T-TO)*(l . + 1 .5*AJ2*(RE/PM)**2* 
*SQRT(OME2)*(l.-3.*BB)) 

AK2= 1.5*AJ2*(RE/PM)**2 
AK32 = .5*AJ3/AJ2 
REO A = RE/AMI 
REOP = RE/PM 
AK42 = AJ4/(AJ2* AJ2) 

S2I = SIN(2.*IMI) 

C2 WS = COS(2 . * APST AR) 


C2AP = COS(2 . * APMI) 

51 = SIN(IMI) 

CI = COS(IMI) 

C2I = COS(2. *IMI) 

CTI = CI/SI 

52 WS = SIN(2 . * APSTAR) 

S2AP = SIN (2 . * APMI) 

SWS = SIN(APSTAR) 

SAP = SIN(APMI) 

CWS = COS( APSTAR) 

CAP = COS(APMI) 

IF(JCIRC.EQ.l) GO TO 1 

ILP=AK2*S2I/(2.-5.*BB)*EMI*EMI*(-7./48. + 5./16.*BB + 35./18.*AK42* 
*(9./56.-3./8.*BB))*(C2WS-C2AP)-AK32*REOP*EMI*CI*(SWS-SAP) 

ELP = -AK2*2.*OME2*SI*SI/(2.-5.*BB)*EMI*(-7./48. + 5./16.*BB + 35./18.* 
*AK42*(9./56.-3./8.*BB))*(C2WS-C2AP) + AK32*REOP*(SWS-SAP)*SI 
RANLP = AK2*.5*CI/(2.-5.*BB)*(EMI*EMI*(-7./12. + 5./2.*BB) + EMI*EMI* 
*BB/(2.-5.*BB)*(-35./12.+25./4.*BB) + 35./18.*AK42*(EMI*EMI*(9./14.- 
*3.*BB) + EMI*EMI*BB/(2.-5.*BB)*(45./14.-15./2.*BB)))*(S2WS-S2AP) + 
*AK32*REOP*EMI*CI/SI*(CWS-CAP) 


I 
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APLP=.5*AK2/(2.-5.*BB)*(BB*(-7./6. + 5./2.*BB) + EMI*EMI*(7./12.-79./ 
*12.*BB + 45./4.*BB*BB) + 4.*BB*(13.-30.*BB)/(2.-5.*BB)*EMI*EMI*(7./48. 
*-5./16.*BB) + 35./18.*AK42*(BB*(9./7.-3.*BB) + EMI*EMI*(-9./14. + 105./ 
*14.*BB-27./2.*BB*BB) + 4.*BB*(13.-30.*BB)/(2.-5.*BB)*EMI*EMI*(-9./56 
*. + 3./8.*BB)))*(S2WS-S2AP) + AK32*REOP*(SI/EMI-EMI*CI*CI/SI)’ l ‘ 
*(CWS-CAP) 

MALP = .5*AK2*SQRT(OME2)/(2.-5.*BB)*(BB*(7./6.-5./2.*BB)*(l.-5./2.*E 
*MI*EMI)-35./18.*AK42*(BB*(9./7.-3.*BB)*OME2))*(S2WS-S2AP)-AK32*REO 
*P*SQRT(OME2)*SI*(L/EMI-EMI)*(CWS-CAP) 

RP = AMF(l.-EMI) 

FE = (RE -B)/RE 

Q0 = RE*FE*SIN(IMI)**2*SIN(APMI)**2/H 
G = DRG*RHONUL*EXP(-(AMI-RP)/H + Q0) 

DB = AMPEMI/H 
DBSQR = DB**2 
DBCUBE = DBSQR*DB 
ESQR = EMI**2 

DD = EL/MMOT*SQRT(OME2)*COS(IMI) 

0 = (RE-B)/H*SIN(IMI)**2 

BO = MBFO(DB) 

B1 = MBFl(DB) 

BB = .5*SIN(IMI)**2 

COMPUTATION OF LONG PERIOD ORBITAL ELEMENTS WITH DRAG 


TDG = TPI/MMOT*( 1 .-AK2*SQRT(OME2)*( 1 .-3. *BB)> 

TPIOT = TPI/(TDG) 

ELM = EL/MMOT 
EP = 1.5*AJ2 

EB = EMI/DB 

DB1 = l./DB 

EBSQR = EB**2 

DEFINITION OF P2K(K = A,E,I,W,M) AND P3I 

Y1A = EB*(4.*(l.-DD)-3.*EB*(17./2.-5.*DD)) 

Z1A = DBl*((l.-DD)**2 + ESQR*(17./2.-5*DD)-8.*EB*(l.-DD) 

1 + 6.*EBSQR*(17./2.-5.*DD)) 

Y1E = DB1*((1.-DD)**2 + ESQR*(1 l./2.-3.*DD)-3.*EB*(l.-DD) 

1 *(3. + .5*DD) 4- 12.*EBSQR*( 1 l./2.-3.*DD)) 

Z1E = DB1*(EMI*(1.-DD)*(3. + DD)-DB1*(2.*(1.-DD)**2 + 5.*ESQR* 
1 (1 1 ./2.-3.*DD)) + 3.*(EMI/DBSQR)*(1.-DD)*(6. + DD)- 
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2 24.*ESQR/DBCUBE*(ll./2.-3.*DD)) 

Y1I = EB*(2.*DD+3.*EB*(.5+1.5*DD)) 

Z1I = DB1*((1.-DD)-ESQR*(.5 + 1.5*DD)-4.*DD*EB-3.*EBSQR* 

1 (l.+3.*DD)> 

Y1ISTR = DBl*(2.*EMI*DD-6./DB*((l -DD)-.75*ESQR*(5.-DD)) + 

1 48.*EMI/DBSQR) 

Z1ISTR = DB1*(U-DD)-.5*ESQR*(11+DD)-4.*EB*(3.+DD) + 

1 12./DBSQR*((1.-DD)-.75*ESQR*(5.-DD))-96.*EMI/DBCUBE) 

Q1A = OME2*(Y1A*BO + Z1A*B1) 

Q1E = OME2*(Y1E*BO + Z1E*B1) 

Q1I = OME2*(Y1I*BO + Z1I*B1) 

Q1ISTR = OME2* ( Y 1 ISTR*BO + Z 1 ISTR*B 1 ) 

= .5*PA-Q1A 
= .5*PE-Q1E 
= . 5* (DPI-PISTAR-2 . *Q 1 1) 

= . 5* ( . 5*PISTAR-Q 1 ISTR + 2 . * Q3I) 

= DBl*((l.-DD)**2/EMI + EMI*(5./2.-DD+ 1.5*DD**2)-3./DB* 

( 1 .-DD)*<3 . + . 5*DD) + 12 . *EMI/DBSQR*(1 l./2.-3.*DD + 
.5*DD**2)) 

= DB1*((1.-DD)*(2. + 1.5 !,! DD)-DB1*(2.*(1.-DD)**2/EMI + EMF 
(43./2.-1 l.*DD + 9./2.*DD**2)) + 6./DBSQR*(l.-DD)* 

(3. + .5*DD)-24.*EMI/DBCUBE*(ll./2.-3.*DD + .5*DD**2)) 

= DBl*((l.-DD)**2/EMI+EMI*(5./2. + .5*DD**2)-3./DB*(l.-DD) 
*(3. + .5*DD) + 12.*EMI/DBSQR*(ll./2.-3.*DD + .5*DD**2)) 

= DB1*((1.-DD)*(2. + 1.5*DD)-DB1*(2.*(1.-DD)**2/EMI + EMI* 
(43./2.-9.*DD + 5./2.*DD**2)) + 6./DBSQR*(l .-DD)* 

(3. + .5*DD)-24.*EMI/DBCUBE*(1 1 ./2.-3.*DD + .5*DD**2)) 

= SQRT(OME2) * ( Y3 W* BO + Z3 W*B 1 ) 

= SQRT(OME2)*(Y3M*BO + Z3M*B1) 

= -Q3W 
= -Q3M 

DEFINITION OF P5K ( K = A,E,I,I*) 

Y5A = 2.0*EMI*(l.-DD**2)-ESQR/DB*(9./2.-DD) 

Z5A = (l.-DD)**2 + ESQR*(1.5 + DD)-EB*(1.-DD)*(3. + DD) + 

1 2.*ESQR/DBSQR*(9./2.-DD) 

Y5E = (l.-DD)**2 + ESQR*(.5-DD)-EB*(l.-DD)*(2. + 1.5*DD) + 

1 3.*ESQR/(2.*DBSQR)*(5.-DD-3.*DD**2) 


P2A 

P2E 

P2I 

P3I 

Y3W 

1 

2 

Z3W 

1 

2 

Y3M 

1 

Z3M 

1 

2 

Q3W 

Q3M 

P2W 

P2M 
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Z5E 


= EMI*(1.-DD)*(1.+3.*DD)-DB1*(1.-DD)*(1.-DD + ESQR* 

1 (3. + 2.*DD)) + EMI/DBSQR*(l.-DD)*(4.+3.*DD)- 

2 3 . * ESQR/DBCUBE* (5 . -DD-3 . * DD* *2) 

Y5I = -2.*EMI*(1.-2.*DD) + ESQR/(2.*DB)*(1. + 7.*DD) 

Z5I = (1 .-DD) + ESQR/2. *(1 .-13. *DD) + EB*(1 ,-3.*DD)- 
1 ESQR/DBSQR*(1 . + 7.*DD) 

Y5ISTR = -2.*EMI*(l.-2.*DD)-DBl*(2.*(l.-DD)-ESQR/2.*(7. + DD)) 
1 + 6. *EMI/DBSQR* ( 1 . + DD) 

Z5ISTR = ( 1 . -DD) + ESQR/2. * (3.-7. *DD) + EB* ( 1 . -7*DD) + 2. /DBSQR 

1 *(2.*(l.-DD)-ESQR/2.*(7. + DD))-12.*EMI/DBCUBE 

2 *(1. +DD) 

P5A = Y5A*B0 + Z5A*B1 

P5E = Y5E*B0 + Z5E*B1 

P5I = Y5I*B0 + Z5I*B1 

P5ISTR = Y 5ISTR* BO + Z5ISTR* B 1 

DEFINITION OF P6K (K = I,W,M) 

Y6I = DB1*((1.-DD)-ESQR/2.*(1. +3.*DD)-3.*EB*(1. +DD) + 

1 6.*ESOR/DBSQR*(l.-DD)) 

U6I = OME2*Y6I 

Z6I = DB1 *(2.*EMI*DD-DB1 *(2.*(1.-DD) + ESQR/2. *(1.-9. *DD)) + 

1 6. *EMI/DBSQR*(1 . + DD)- 12.* ESQR/DBCUBE* ( 1 . -DD)) 


V6I = OME2*Z6I 

Y6W = DB1*((1.-DD)*(2. + 1.5*DD)-3.*EB*(5./2.-.5*DD)) 

U6W = SQRT(OME2)*Y6W 

Z6W = DBl*((l.-DD)**2/EMI + EMI*(5./2.-DD+ 0.5*DD**2)- 
1 2./DB*(l.-DD)*(2. + 1.5*DD) + 6.*EMI/DBSQR*(5./2.-.5*DD)) 

V6W = SQRT(OME2)*Z6W 

U6M = U6W 

Z6M = DBl*((l.-DD)**2/EMI + EMI*(5./2.-.5*DD**2)-2./DB*(l.-DD) 

1 * (2. + 1.5*DD) + 6.*E/DBSQR*(5./2.-.5*DD)) 

V6M = SQRT (OME2) * Z6M 
P6I = U6I*B0 + V6I*B1 
P6W = U6W*B0 + V6W*B1 
P6M = U6M*B0+V6M*B1 
C 

C DEFINITION OF AWBAR,A\A0\B\B0\A",A0",AC\AC0\BS\BS0’ 

C 

C ALDOT = MEAN RATE OF THE SUN RELATIVE TO THE EARTH 
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ALI = INITIAL MEAN LONGITUDE OF THE SUN AT TIME TO 
EPS1 = THE OBLIQUITY OF THE ECLIPTIC 
DLMDA = ANGLE IN RIGHT ASCENSION THAT CENTER 
OF BULGE LAGS THE SUN 
EPS1 = 23.44D0*RPD 
ALDOT = (0.9856D0*RPD)/86400.D0 
DLMDA = 30.D0*RPD 
SEPS1 = SIN(EPSl) 

COISQ = COS(IMI/2.)**2 
SISQ = SIN(IMI/2.)**2 
APDT = AK2*MMOT*(2.-5./2.*SIN(IMI)**2) 

RADT = -AK2*MMOT*COS(IMI) 

ELI - ALDOT + APDT 
EL2 = ALDOT - APDT 
EL3 = ALDOT-RADT + APDT 
EL4 = ALDOT-RADT-APDT 
EL5 = ALDOT + 3.* APDT 
EL6 = ALDOT - 3.*APDT 
EL7 = ALDOT-RADT +3.* APDT 
EL8 = ALDOT -RADT -3 . * APDT 
SLW = ALI + APMI + ELI *(T-TO) 

SLWO = ALI + APMI 

DLW = ALI-APMI + EL2*(T-TO) 

DLWO = ALI-APMI 

SLOW = ALI-RANMI + DLMDA + APMI + EL3 * (T-TO) 

SLOWO = ALI-RANMI + DLMDA + APMI 

DLOW = ALI-RANMI + DLMDA-APMI + EL4* (T-TO) 

DLOWO = ALI-RANMI + DLMDA- APMI 
SL3W = ALI+3.*APMI + EL5*(T-TO) 

SL3W0 = ALI + 3.*APMI 

DL3W = ALI -3 . * APMI + EL6* (T-TO) 

DL3W0 = ALI-3.*APMI 

SL03W = ALI-RANMI + DLMDA +3.* APMI + EL7* (T-TO) 

SLO3W0 = ALI-RANMI + DLMDA + 3.* APMI 
DL03W = ALI-RANMI + DLMDA-3. *APMI + EL8*(T-TO) 

DLO3W0 = ALI-RANMI + DLMDA-3 . * APMI 


SNA 

= SIN(SLW) 

SNAO 

= SIN(SLWO) 

SNB 

= SIN(DLW) 

SNBO 

= SIN (DLWO) 

SNC 

= SIN(SLOW) - 

SNCO 

= SIN(SLOWO) 

SND 

= SIN(DLOW) 
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SNDO = SIN(DLOWO) 
SNE = SIN(DL3W) 
SNEO = SIN(DL3W0) 


SNF = SIN(SL3W) 

SNFO = SIN(SL3W0) 

SNG = SIN(DL03W) 

SNGO = SIN(DLO3W0) 

SNH = SIN(SL03W) 

SNHO = SIN(SLO3W0) 

CSA = COS(SLW) 

CSAO = COS(SLWO) 

CSB = COS(DLW) 

CSBO = COS(DLWO) 

CSC = COS(SLOW) 

CSCO = COS(SLOWO) 

CSD = COS(DLOW) 

CSDO = COS(DLOWO) 

AWBAR = (RE/PM)**2*MMOT*(2.-5./2.*SIN(IMI)**2) 
APRIM = .5*SEPSI*SIN(IMI)*(EP/EL2*SNB-EP/EL1*SNA) 

1 + SISQ*EP/EL3*SNC + COISQ*EP/EL4*SND 

BPRIM = -.5*SEPSI*SIN(IMI)*(EP/EL2*CSB-EP/EL1*CSA) 

1 + SISQ*EP/EL3*CSC-COISQ*EP/EL4*CSD 

ACPRM = EP/4.*SEPSI*SIN(IMI)*(1./EL6*SNE-1./EL5 

1 *SNF-1./EL2*SNB+1./EL1*SNA) + 

2 EP/2.*COISQ*(l./EL8*SNG + l./EL3*SNC) 

3 + EP/2.*SISQ*(1./EL7*SNH + 1./EL4*SND) 

BSPRM = EP/4.*SEPSI*SIN(IMI)*(1./EL6*SNE-1./EL5 

1 *SNF + 1 ./EL2*SNB-1 ./ELI *SNA) + 

2 EP/2.*COISQ*( 1 ./EL8*SNG- 1 ./EL3*SNC) 

3 + EP/2.*SISQ*(1./EL7*SNH-1./EL4*SND) 

ADLPM = -EP**2/2.*SEPSI*SIN(IMI)*(1./EL2**2*CSB-1./EL1 

1 **2*CSA)-SISQ*EP**2/EL3**2*CSC- 

2 COISQ*EP**2/EL4**2*CSD 

AOPRIM = .5*SEPSI*SIN(IMI)*(EP/EL2*SNB0-EP/EL1*SNA0) 
1 + SISQ*EP/EL3*SNC0 + COISQ*EP/EL4*SNDO 

BOPRIM = -.5*SEPSI*SIN(IMI)*(EP/EL2*CSB0-EP/EL1*CSA0) 
1 + SISQ*EP/EL3*CSCO-COISQ*EP/EL4*CSDO 

ACOPRM = EP/4.*SEPSI*SIN(IMI)*( 1 ./EL6*SNE0-1 ./EL5 
1 *SNF0-1 ./EL2*SNB0 + 1 ./EL1*SNA0) + 
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2 EP/2.*COISQ*(1./EL8*SNGO+ l./EL3*SNC0) 

3 + EP/2.*SISQ*(1./EL7*SNH0 + l./EL4*SND0) 

BSOPRM = EP/4.*SEPSI*SIN(IMI)*(1./EL6*SNE0-1./EL5 

1 *SNFO + 1 . /EL2*SNB0- 1 . /EL 1 *SNAO) + 

2 EP/2. *COISQ*(l ,/EL8*SNG0- 1 ./EL3*SNC0) 

3 + EP/2. *SISQ*(1 ./EL7*SNHO-l ./EL4*SND0) 

AODLPM = -EP**2/2.*SEPSI*SIN(IMI)*(1./EL2**2*CSB0-1./EL1 

1 **2*CSA0)-SISQ*EP**2/EL3**2*CSC0- 

2 COISQ*EP**2/EL4**2*CSDO 

DAPM = APMI + APDT* (T-TO) 

DS2W = SIN(2.*DAPM) 

DS2WI = SIN(2.*APMI) 

DS4W = SIN(4.*DAPM) 

DS4WI = SIN(4.*APMI) 

DC2W = COS(2.*DAPM) 

DC2WI = COS(2.*APMI) 

DSI = SIN(IMI) 

DSISQ = DSI*DSI 
GOEP = G/EP 
REOPSQ = (RE/PM)**2 

ALPD = -TPIOT* (2 . * GOEP* AMI **2*(Q*P2A/(2.* A WB AR) * 

1 (DS2W-DS2WD+ F*P5A*(APRIM-A0PRIM))) 

ELPD = -TPIOT*(2.*GOEP*AMI*OME2*(Q*P2E/(2.*AWBAR)* 
1 (DS2W-DS2WI) + F*P5E*(APRIM-A0PRIM))) 


DILP = -TPIOT*(.5*GOEP*AMI*ELM*DSI/SQRT(OME2)*((PISTAR + Q* 

1 P2I)/(2.*AWBAR)*(DS2W-DS2WI) + Q*P3I/(4.*AWBAR)* 

2 (DS4W-DS4WI) + F* (P5I* (APRIM-AOPRIM) + P5ISTR* 

3 (ACPRM-AC0PRM)-2.*P6I*(BSPRM-BS0PRM)))) 

RANLPD = TPIOT*(GOEP*REOPSQ*MMOT*AMI*COS(IMI)*(Q/(4. 

1 *AWBAR**2)*(7.*P2A-8.*EMI*P2E)*(DC2W-DC2WI)- 

2 F*(ADLPM-A0DLPM)*(7.*P5A-8.*EMI*P5E))) 

APLPD = TPIOT*(GOEP*AMI*SQRT(OME2)*(Q*P2W/AWBAR* 

1 (DC2W-DC2WI)-2.*F*P6W !|! (BPRIM-BOPRIM)))-TPIOT* 

2 (GOEP*REOPSQ*MMOT*AMI*(2.-5./2.*DSISQ)*(Q/(4.* 

3 AWBAR**2)*(7.*P2A-8.*EMI*P2E) :,: (DC2W-DC2WI)- 

4 F*(ADLPM-A0DLPM) :t: (7.*P5A-8.*EMPP5E))) 

DMALP = -TPIOT*(GOEP*AMI*(Q*P2M/AWBAR*(DC2W-DC2WI)- 
1 2.*F*P6M*(BPRIM-BOPRIM)))-TPIOT*(GOEP*REOPSQ* 
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2 MMOT*AMI*SQRT(OME2)*(l.-1.5*DSISQ)*(Q/(4.* 

3 AWBAR**2)*(7.*P2A-6.*EMI*P2E)*(DC2W-DC2WI)-F* 

4 (ADLPM-A0DLPM)*(7.*P5A-6.*EMI*P5E))) 

ILP = ILP + DILP 

ELP = ELP + ELPD 
APLP = APLP + APLPD 
RANLP = RANLP + RANLPD 
MALP = MALP + DMALP 
RETURN 

1 ETAB = QOB*ALP*SIN(BO*(T-TO) + ABET) + EOB 
ETAI = QOB*ALP*SIN(ABET) + EOB 
ZMB = ALP* COS(BO* (T-TO) + ABET) 

RANLP = .5*(AJ3/AJ2)*REOA*CTI*(ZMB-ALP*COS(ABET)) 

LAMLP = -.5*(AJ3/AJ2)*REOA*C2I/SI*(ZMB-ALP*COS(ABET)) 

ILP = -.5*(AJ3/AJ2)*REOA*CI*(ETAB-ETAI) 

RETURN 

END 
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